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On  singular  values  of  Hankel  operators  of  finite  rank 

W.  B.  Graggf  and  L.  Reichel  % 


Abstract 

Let  H  be  a  Hankel  operator  defined  by  its  symbol  p  =  x/x  where  x  is  a  monic  polynomial  of  degree 
n  and  x  is  a  polynomial  of  degree  less  than  n.  Then  H  has  rank  n.  We  derive  a  generalized  Takagi 
singular  value  problem  defined  by  two  n  x  n  matrices,  such  that  its  n  generalized  Takagi  singular 
values  are  the  positive  singular  values  of  H.  If  p  is  real,  then  the  generalized  Takagi  singular  value 
problem  reduces  to  a  generalized  symmetric  eigenvalue  problem.  The  computations  can  be  carried 
out  so  that  the  Lanczos  method  applied  to  the  latter  problem  requires  only  0(n  log  n)  arithmetic 
operations  for  each  iteration.  If  x  and  x  ,drt  given  in  power  form,  then  the  elements  of  all  n  x  n 
matrices  required  can  be  determined  in  0(n7)  arithmetic  operations. 
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1.  Introduction 

Let  H  —  [f?y+k]^fc=o  be  a  be  a  Hankel  operator  defined  by  its  rational  symbol  p  =  n/x,  where 


n-l 


tt(A)  :=  J2  ***'     and     X(A)  :=  ]£»*'.  Xn  =  1.  (1.1) 

J=0  j=0 


We  assume  that  n  and  x  have  no  common  zeros.  The  elements  rjj  of  H  are  then  given  by 


In  order  to  simplify  our  presentation,  we  assume  that  the  zeros  {Afc}£=1  of  x  are  distinct.  How  our 
formulas  need  to  be  modified  in  order  to  remove  this  assumption  is  discussed  in  Remark  1.1  below. 
Hence  p  has  a  partial  fraction  decomposition 

k=i  * 

Expansion  of  the  right  hand  side  of  (1.3)  into  a  geometric  series,  and  comparison  with  (1.2),  yields 

n 

!?y  =  £c*Aj.  (1-4) 

We  now  express  (1.4)  in  matrix  form.  Let 

A:  =  diag[ai,a2,...an]  eCnXn,  (1.5) 

A:=   dioy[A1,A2,...Ari]e<rnXnJ  (1.6) 

and  introduce  the  Vandermonde  matrix 

Vb  :  =  [Ai+JJjio  .«  €7»x».  (1.7) 

Define 

^  :=  [VylJLo  c  CcoXn,  (1-8) 

where 

Vy:=VbA^,  j>l.  (1.9) 

Then  (1.4)  can  be  written  as 

i/  =  VA7T.  (1.10) 

Let  I2  denote  the  vector  space  <C°°  equipped  with  the  Euclidean  norm. 


Proposition  1.1.  H  :  I2  ->  I2  bounded  <&  \Xk\  <    1    for    1  <  k  <  n. 

Proof.    The  proposition  holds  independent  of  the  multiplicity  of  the  A^.    In  the  present  proof  we 
assume  that  the  A^  are  distinct.  The  proof  for  confluent  A*,  is  commented  on  in  Remark  1.1. 

Let  ei  =  [ej]?l0  c  <C°°  be  the  axis  vector  with  e0  =  1-  Then 

h  =  [»7y]yl0  :  =  Hex  €.  I2  =>  rj3  ->  0  as  j  ->  oo  => 

\Xk\    <    1  for  1  <  k  <  n, 
where  the  last  implication  follows  from  (1.4). 

Conversely,  assume  that  |Afc|  <    1  for  1  <  k  <  n.  Then  by  (1.8)  -  (1.10)  we  obtain 


ll^ll2<||A|!2|^||2<||A||2||V0||l||^(A^Ar||? 

=  l|A||2||^!l^||(/-(A//A)ri)-1||l        - 
We  assume  henceforth  that  \\k\    <    1  for  1  <  k  <  n.  Introduce 

U:=VV-0\  (1.11) 

Ho-VqAV*.  (1.12) 

Then  Ho  has  rank  n.  We  note  ,  by  comparing  (1.12)  with  (1.10),  that  Hq  is  the  leading  principal 
n  x  n  submatrix  of  //.  From  (1.10)  -  (112)  it  follows  that 

H  =  UH0UT.  (1.13) 

The  leading  n  x  n  submatrix  of  U  is  In,  the  n  x  n  identity  matrix.  U  therefore  is  of  rank  n  and 
can  be  factored 

U  =  QR,     Qe(U°°Xn,      Re(CnXn, 

where  QH Q  —  In  and  R  is  a  nonsingular  right  triangular  matrix.   We  obtain 

a+{H)  =  a+{QRH0RTQT)  =  a{RH0RT),  (1.14) 

where  a  denotes  the  set  of  singular  values  and  a+  denotes  the  subset  of  the  positive  ones. 

The  n  x  n  matrix  RH0RT  is  complex  symmetric.  Takagi  [Tal],  [Ta2]  showed  the  existence  of  a 
complex  symmetric  singular  value  decomposition 

RH0RT  =  WEWt,  We(D,,Xn,  X  =  diag[aJ}(T2,-  <?n},  (1.15) 

where  W^W  =  /„  and  a }  >  0  are  the  singular  values  of  RHoRT .  In  Section  2  we  present  an 
elementary  proof  of  the  existence  of  this  decomposition.  Let  W  —  [wi,W2,  ■■■  u^ri ] ,  w}  (.  Cn .  Then 
(1.15)  can  be  written  as  the  Takagi  singular  value  problem 

RHQRTw3  =  Wjaj,      w'/wk  =  5:h,      1  <  j,k  <  n,  (1.16) 
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where  the  bar  denotes  complex  conjugation  and  80k  is  Kronecker's  8  function.  The  problems  (1.15) 
-  (1.16)  could  be  solved  by  the  algorithm  described  in  [BGG],  but  this  would  require  RH0RT  to  be 
explicitly  computed.  In  order  to  avoid  these  matrix  multiplications  we  let  v3  :=  RH w3  and  obtain 
from  (116)  the  generalized  Takagi  singular  value  problem 

H0vj  =  {RHR)-1v,a,,     vf{RHR)~1vk  =  8jk,  l<j,k<n.  (1.17) 

The  solution  of  (1.17)  requires  (RH  R)'1  to  be  known.  In  Section  3  we  show  that 

{R'JR)-1  =  I-B0B»,  (1.18) 

where  Bo  c  CnXn  is  a  triangular  Toeplitz  matrix.  The  elements  of  Bo  and  Ho  can  be  determined 
from  the  coefficients  of  n  and  x  in  0(n  log  n)  arithmetic  operations  by  the  fast  Fourier  transform 
(FFT)  method.  This  is  demonstrated  in  Section  4.  Section  5  shows  that 


RHR  =  T1M0T{I}     TuM0eCnXn,  (1.19) 

where  T\  and  Mo  are  Toeplitz  matrices,  and  describes  a  numerical  scheme  for  the  computation 
of  this  factorization  from  (1.16)  in  0(n2)  arithmetic  operations.  We  also  present  a  Hermitian 
factorization  of  RH  R  into  n  x  n  triangular  matrices. 

The  factorization  (1.19)  may  be  of  interest  for  the  numerical  solution  of  (1.17).  Assume  that 
the  coefficients  of  n  and  x  are  real  valued.  Then  Ho,  (RH R)~1  e  IRnXn,  and  (1.17)  reduces  to 
a  generalized  symmetric  eigenvalue  problem.  The  Lanczos  method  ([Pa,  Section  15.11],  [ER]) 
would  appear  suitable  for  solving  this  eigenproblem  for  the  following  reason.  Let  C  e  <CnXn  be  a 
Hankel  or  Toeplitz  matrix  and  let  v  e  <Vn  be  arbitrary.  It  is  well  known  that  Cv  can  be  computed 
in  0{n  log  n)  arithmetic  operations  using  FFTs.  Hence  H0v,  (RH R)~1v  and  (RH R)v  can  be 
computed  in  0(n  log  n)  arithmetic  operations,  where  we  use  (118)  -  (1.19).  Each  iteration  of  the 
Lanczos  algorithm  given  in  [Pa,  p. 324]  therefore  requires  only  0(n  log  n)  arithmetic  operations. 

The  computation  of  singular  values  of  H  is  important  in  Hankel  norm  approximation  problems  of 
systems  theory,  such  as  the  model  reduction  problem  [Gil.  The  approximation  of  functions  by  the 
Caratheodory  -  Fejer  method  yields  another  application  [Gu],  [Tr]. 

Other  methods  for  reducing  the  singular  value  problem  for  H  to  a  finite  dimensional  one  have  been 
described  by  Kung  and  Gutknecht  [Gu]  and  Young  [Yo].  These  methods,  however,  do  not  preserve 
symmetry.  Moreover,  Young's  approach  requires  generally  0(n3)  arithmetic  operations  to  compute 
the  matrices  required. 

Remark  1.1.  Formulas  (1.3)  -  (1.8)  and  the  proof  of  Proposition  1.1  require  distinct  A^.  This 
restriction  can  be  removed.  Assume  first  that  Aj  =  A2  =  ...  =  An.  Then  (1.3)  -  (1.4)  have  to  be 
replaced  by 

n  oo      .  , 

k=l  j  =  o 

In  (1.5)  A  has  to  be  substituted  by  the  upper  triangular  Hankel  matrix 

i4=[ai+Jfc+1£fc*0e  CnXn;      ap  :=  0,      p  >  n. 
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The  matrix  A  in  (16)  has  to  be  replaced  by  the  Jordan  matrix  with  all  diagonal  elements  equal  to 
A:  and  all  superdiagonal  elements  equal  to  one.    The  matrix  V0  in  (1.7)  need  be  replaced  by  the 
onfluent  Vandermonde  matrix.  For  instance,  we  obtain  for  n  =  3 


A  = 


A!         1 

Ai       1 


Vn  = 


Ai        1 
A?     2A, 


With  A,  A  and  V0  modified  as  described,  we  define  V:  and  V  by  (18)  -  (1.9),  U  by  (111)  and 
Hq  by  (1.12).  Then  (1.10)  and  (1-13)  hold  and  Hq  is  the  leading  principal  n  x  n  submatrix  of  H . 
Also  (1.14)  -  (119)  remain  valid.  Proposition  1.1  can  be  shown  by  replacing  (1.4)  by  (14'),  and 
by  bounding  the  sum 


nE(A"A) 


«>wi 


j=0 


where  A  now  is  a  Jordan  matrix.    This  sum  is  bounded  if  |A] 
valid. 


<  1,  and  the  proposition  remains 


In  general,  when  the  A*,  are  of  arbitrary  multiplicity,  A  in  (1.5)  has  to  be  replaced  by  a  block 
diagonal  matrix,  where  each  block  is  an  upper  triangular  Hankel  matrix.  The  blocks  are  of  the 
same  sizes  as  the  multiplicities  of  the  A^,  and  the  number  of  blocks  equals  the  number  of  distinct  A^ . 
A  in  (1.6)  is  replaced  by  a  Jordan  matrix  with  Jordan  boxes  of  the  same  sizes  as  the  multiplicities 
of  the  Afc,  and  the  number  of  boxes  equal  to  the  number  of  distinct  A^.  Vq  in  (1.7)  is  replaced  by 
an  appropriate  confluent  Vandermonde  matrix.  With  these  changes  (1.10)  -  (1.19)  are  valid,  and 
so  is  Proposition  1.1.  We  omit  the  details  since  the  numerical  computations  are  independent  of  the 
multiplicity  of  the  A*.,     m 


2.  The  Symmetric  Singular  Value  Decomposition 


In  this  section  we  present  an  elementary  proof  of  Takagi's  theorem,  i.e.  we  show  the  existence  of 
a  symmetric  singular  value  decomposition  of  a  complex  symmetric  matrix.  Let  C  —  CT  e  CnXn , 
and  define  A,  B  e  RnXn  by  C  :=  A  +  iB,  i  :=  \/^\.  Then  A  =  ATand  B  =  BT ,  so  the  matrix 


C  := 


A        B 
B     -A 


is  real  and  symmetric.  Let  {0]}Tj=\  be  the  positive  eigenvalues  of  C  and  form 

£  :=  diag{aua2,...,ar}. 
Let 

with 

and 

Write  (2.1)  as 

and  note  that  (2.2)  also  can  be  written  as 

AV     + 

BV  U){-E), 

(-S) 


A        B 

U 

U 

■s 

B     -A 

V 

V 

2_ 

U,Ve  Rnxr 

UTU  +  VTV  =  Ir. 

AU     +     BV            UE 

.BU     - 

AV 

VH 

(AV     +     B{-U)  V{-X) 


i.e. 


A 
B 


B 
A 


V 

-u 


V 

u 


with 


VTV  +  {-U)T{-U)  =  Ir. 


(2.1) 


(2.2) 


(2.3) 


Hence  C  has  at  least  r  negative  eigenvalues.  We  could  also  have  let  ay  be  the  negative  eigenvalues  of 
C  and  then  (2.3)  would  have  given  us  positive  ones.  We  therefore  may  assume  that  ±aj ,  ±02,  •••>  i^Y 
are  all  the  nonzero  eigenvalues  of  C . 

Since  eigenvectors  associated  with  distinct  eigenvalues  of  a  real  symmetric  matrix  are  orthogonal, 
we  have 


0  =  [vT,-u: 


=  VTU  -  UTV. 


The  spectral  resolution  of  C  is  thus 


A         B 
B      -A 


U 
V 


V 

u 


UT 
VT 


VT 
UT 


which  yields 

f  A  =  UEUT  -VEVT 
\B  =  VEUt  +  UEVt. 

Therefore 

C  =  A  +  iB  =  UZUT  -  VEVT  +  i{VT,UT  +  UY,VT) 


=  {U  +  iV)E(£/T  +  tVr)  =  WEWT  =  ^  akww 
U  +  iV  =:W  =  {wuw2,...,wr},  wk(  <Dn 


T 

h  > 


fc=l 

where 


Moreover 

WHW  =  (UT  -  iVT)(U  +  iV)  =  (UTU  +  VTV)  +  i(UTV  -  VTU)  =  Ir. 

If  r  <  n  then  one  may  replace  £  by 

S0  ~diag\a1,a2l...,aT,0>...,0}  c  RnXn 

and  W  by 

W0  :=  [wi,w2,  ...,wr,wr+1,  ...wn]  e  €nXn, 

where  wr+i,  ...,wne  Cn  are  chosen  so  that  WqWq  =  In.    m 


3.  A  Simple  Expression  for  (RH R)    * 

In  this  section  we  derive  (1.18).  Introduce  the  Frobenius  matrix 

F:={e2,e3}...,en>-f}eCnXn, 

where 

ey:=[*iy,fci,...,$Bjfretfn.    2  <  3  <  ",  (3-1) 

f:  =  {Xo,Xu-,Xn-i]TcCn- 

Then  F  is  the  companion  matrix  of  x  and 

FTVQ  =  V0A.  (3.2) 

Throughout  this  section  V0  and  A  are  defined  by  (1.6)  -  (1.7)  if  the  A*,  are  distinct.  For  confluent 
Afc  we  modify  V0  and  A  according  to  Remark  1.1.  The  following  lemma  shows  that 


G  :=  RHR  (3.3) 

satisfies  a  Stein  equation.  This  will  enable  us  to  obtain  a  simple  expression  for  G_1  by  an  application 
of  the  Sherman-Morrison- Woodbury  formula. 

Lemma  3.1.  G  is  the  unique  solution  of  the  Stein  equation 

X  -  FnXFnH  =  In,  Xe  Cnxn.  (3.4) 

Proof.  By  (1.8),  (1.9)  and  (1.11)  we  obtain 

oo 

RHR  =  UHU  =  J2Vo~H{knk)HvoHvoknkVo~\  (3-5) 

and  (3.2)  yields  now 

CO 

G  =  J2Fnk{Fnk)H.  (3.6) 

k  =  Q 

The  series  in  (3.5)  -  (3.6)  converge  because  |Afc|  <  1  for  all  k.  Substitution  of  (3.6)  into  (3.4) 
shows  that  G  solves  (3.4).  The  unicity  follows  from  |Afc|  <  1  for  all  k.  The  latter  can  be  seen  by  a 
similarity  transform  of  Fn  to  Schur  triangular  form,    m 

Introduce  the  cyclic  downshift  operator  in  <C2n 

JB:=[e2,e3....,cn,Cl]  e  C2nx2n , 

where 

tj  :=l6l3i62j;...,62ntJ]T  eIR2n.  (3.7) 

Let 

t:=[xo,Xi,-,Xn,0,0,...,0)Te€2n, 
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and  define  the  Toeplitz  matrix  T  of  parallelogram  form 

T  :=  \t,Et,E2t,...,En-h}  e  C2nxn. 


(3.8) 


Let  T0  be  the  leading  n  x  n  submatrix  of  T,  and  let  Tj  be  the  trailing  n  x  n  submatrix  of  T.  Then 
T0  is  a  left  triangular  Toeplitz  matrix,  and  Ti  is  a  unit  right  triangular  Toeplitz  matrix. 

Example  3.1.  Let  n  =  3.  Then 

Xo 

Xi     Xo 

X2     X2     Xo 

X3       X2       Xl 

X3       X2 

X3 


To 


Xo 

X3 

X2 

Xi 

X\ 

Xo 

Ti  = 

X3 

X2 

X2 

Xi 

Xo. 

X3_ 

where  we  note  that  X3  —  1-     • 

Lemma  3.2.    Let  To  and  T\  be  defined  as  above.  Then 

T0"  T0  +  T**TX  =  TqTo"  +  T^Ti 


a 


(3.9) 


Proof.    Let  N  :=  THT  =  T^T0  +  T?TX.    We  first  show  that  N  is  a  Toeplitz  matrix.    Let  ey  be 
defined  by  (3.1).  Then  by  (3.8)  we  have  for  1  <  j,k  <  n, 

ej  Nek  =  e0THT  ek  =  tH{EHy-lEk~H  =  tHEk~>t, 
where  we  have  used  that  E     =  E~l .  We  next  define  the  reversal  matrix 

Toeplitz  matrices  are  counter  symmetric,  i.e.  N  =  J NT  J.  Using  that  N  is  counter  symmetric  and 
Hermitian  yields 

T0"To  +  T?Tx  =  N  =  JNTJ  =  JNJ  =  J{Tfn  +  T^TX)J 

=  JT^J  ■  JTqJ  +  JTjj  •  JYXJ  =  T0T^  +  TXT? .      , 

The  next  lemma  presents  a  Gaussian  factorization  of  Fn  in  terms  of  To  and  T\.  This  will  be  used 
together  with  Lemma  3.1  to  express  G~l  in  terms  of  T0  and  Tj . 


Lemma  3.3. 


-  i 


Fn  =  -ToTf1. 


(3.10) 


Proof.  We  first  show  that 


T    rrT] 


Po,^ 


V0 


V0A' 


0. 


Let  ey  be  defined  by  (3.7)  and  assume  for  the  moment  that  the  A^  are  distinct.  Then 

V0 


$\tLiT] 


V„Af 


ek  =  XMK 


j-i 


(3.11) 


(3.12) 


and  the  right  hand  side  vanishes  for  1  <  j,  k  <  n.  If  the  A^  are  confluent,  then  the  right  hand  side 
expression  of  (3.12)  contains  derivatives  of  x(^)  evaluated  at  A*..  The  right  hand  side  of  (3.12), 
however,  still  vanishes  and  (3.11)  holds. 

We  now  write  (3.11)  as 

T0rVo  +  T1TVoAn  =  0 

and  apply  (3.2).  This  shows  (3.10).      • 

We  are  now  in  a  position  to  show  (1.18).  By  (3.4)  G  satisfies 

G  =  I  +  FnG  FnH 
and  an  application  of  the  Sherman-Morrison- Woodbury  formula  yields 

G~l  =  (/  +  FnG  F"")-1  =  /  -  Fn{G~1  +  pnH Fn)~lFnH .  (3.13) 

We  now  determine  an  expression  for 

Y  —  I-G'1.  (3.14) 

Substitute  Y  and  (3.10)  into  (3.13)  to  obtain 

Y  =  To  (To"  To  +  T^Ti  -  T^YT^T^ .  (3.15) 

In  order  to  determine  a  simple  expression  for  Y  from  (3.15)  we  need  the  following  observation, 
which  is  also  central  to  Section  4.  To  and  T£~  are  both  left  triangular  n  x  n  Toeplitz  matrices. 
Multiplication  of  To  with  T^H  can  be  identified  with  polynomial  multiplication,  see  [Hel,  Section 
1.3]  and  Section  4.  Since  multiplication  of  polynomials  commutes,  we  obtain 

T0T{H  =T{HT0.  (3.16) 

From  the  correspondence  between  polynomials  and  left  triangular  Toeplitz  matrices  it  also  follows 
that  TqTi       is  a  left  triangular  Toeplitz  matrix. 

Lemma  3.4.  Equation  (3.15)  has  the  unique  solution 

Y  =  Tr^ToTo^Tj-1  =  ToT^Tr1^".  (3.17) 

Proof.  Unicity  follows  from  (3.14)  and  that  (3.4)  has  a  unique  solution.  From  (3.16)  we  obtain 

T~HT0T^T^  =T0TrHT^1Tf.  (3.18) 

Now  substitute 

Y  ^T^TqT^T-1 

into  (3.15).  We  obtain 

T^HT0TgT^1  =  To(To"To  +  TfT,  -  T0T^)-XT^ .  (3.19) 

10 


An  application  of  (3.9)  reduces  (3.19)  to  (3.18).    The  latter  has  already  been  shown  to  be  valid. 
Therefore  (3.17)  solves  (3.15).     ■ 

Let 

Bo  :=  ToT~T  =  TrTT0.  (3.20) 

Then  Bq  is  a  left  triangular  n  x  n  Toeplitz  matrix.  By  (3.14)  and  (3.17) 

G       =  /  —  BqB0  =  I  —  Bq  B0. 

From  (3.3)  it  now  follows  that 

{RHR)~1  =  I-  B0B^.  (3.21) 
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4.  Computation  of  Hq  and  Bq 

We  summarize  some  results  in  [He  1,  Section  1.3]  and  [He  2,  Section  13.9]  in  order  to  show  that  the 
elements  of  Ho  and  Bo  can  be  computed  in  0(n  log  n)  arithmetic  operations  from  the  coefficients 
Xj  of  x  and  n3  of  7r,  see  (l.l).  To  a  polynomial  or  power  series 

n-l 

we  associate  the  left  triangular  n  x  n  Toeplitz  matrix 

^  =  [&-fc]?,fc=o.         ?y=0fori<0, 

and  we  write  f  — *  Z .  If  f  (A)  is  a  polynomial  and  X  a  left  triangular  n  x  n  Toeplitz  matrix  such 
that  £  — »  X,  then  it  is  easily  seen  that  c£  — »  ZX.  In  particular,  ZA  is  a  left  triangular  n  x  n 
Toeplitz  matrix.  From  £f  =  $£  and  £c  — *  XZ  is  follows  that  ZA"  —  XZ . 

Assume  that  f o  ^  0  and  let  1/c  -»  Z'.  Then  1/f  •  f  -»  /,  Z'Z  and  ZZ\  We  obtain  Z'  =  Z_1  and 
therefore  Z_1  is  a  left  triangular  Toeplitz  matrix. 

Example  4.1.  We  have  x  ~ *  To.  Let 

n 

X(A)  :=  Art  x(l/A)  =  2>B-yA'.  (4.1) 

Then  x  ~ *  3Ti    and  the  Blaschke  product 

4  ^T0T^H  =  B0.  (4.2) 

X  ■ 

Now  let  £(A)  and  c(A)  be  arbitrary  polynomials  such  that  f  (0)  7^  0.  Henrici  [He2,  Theorem  13. 9e  ] 
shows  that  the  first  n  coefficients  in  the  MacLaurin  expansion  of  f(A)/c(A)  can  be  computed  in 
0(n  log  n)  multiplications.  The  proof  uses  FFT.  It  is  easily  seen  that  the  number  of  additions  also 
is  0{n  log  n). 

From  Xn  —  1  and  (4.1)  we  obtain  x(0)  7^  0.  Hence,  the  first  n  terms  in  the  MacLaurin  expansion 
of  x/x  °an  be  computed  in  0(n  log  n)  arithmetic  operations.  By  (4.2)  therefore  T0T^H  =  Bo  can 
be  computed  in  0(n  log  n)  arithmetic  operations. 

Because  A"x(l/A)  7^  0  for  A  =  0,  we  can  compute  the  first  n  terms  in  the  MacLaurin  expansion  of 

A",(l/A)      yyv+.     0(jn 

in  0(n  log  n)  arithmetic  operations.  This  shows  that  i/0  can  be  computed  in  0(n  log  n)  arithmetic 
operations. 
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5.  A  Factorization  of  RH  R 

It  follows  from  (3.3)  and  (3.20)  -  (3.21)  that 


G'1  =  (RHR)      =I-B0  Bg  =  1-  Tf^To^rf1, 


and  therefore 


T^G'1^  =  T^TX  -  T0T**  =:  M0" 


i 


(4.1) 
(4.2) 


The  expression  defining  M0_1  is  a  Gohberg-Semencul  formula  for  the  inverse  of  an  n  x  n  Toeplitz 
matrix,  see,  e.g.,  [Io,  Theorem  18.2,  p.  152].  We  denote  this  Toeplitz  matrix  by  Mo.  From  the  left 
hand  expression  of  (4.2)  and  the  nonsingularity  of  T^  and  R  it  follows  that  M0  is  Hermitian  and 
positive  definite.  The  desired  factorization  of  RH  R  is 


>// 


H 


R"  R  =  T1  M0  T( 

We  will  now  show  how  Mq  can  be  computed.    The  computation  involves  running  the  Levinson 
algorithm  backwards. 

Consider  the  related  Gohberg-Semencul  formula,  see,  e.g.,  [Io,  Theorem  18.1,  p.  148]  or  [AG], 

H 


m^  = 


Xn        X'i  — 1 


Xo 


Xn-l 

Xn 


Xn-l 


0 
XO 

\1 


lXn-1 


0 

Xo 


Xi      Xo     OJ  lXn-1 


Xo 


Xn-l 
Xn 


H 


(4.3) 


Xi      Xo     0 


where  the  four  triangular  Toeplitz  matrices  define  the  inverse  of  an  (n  +  1)  x  (n  +  1)  Hermitian 
Toeplitz  matrix.  Denote  this  Toeplitz  matrix  by  Mx .  Then  M0  is  the  leading  principal  n  x  n 
submatrix  of  M\,  see  [Io,  Theorems  18.1  -  18.2]. 


Let  Rx  :=  [pjk]^k=o     e  <T(n  +  1)x(n+1)  be  the  unit  right  triangular  matrix,  and  let 
Di  :=  diag\8(),5it  ...,8n]  be  the  diagonal  matrix  such  that 

R?  Mj  R1  =  Dl. 


(4.4) 


L3 


Given  Mi  =  [/iy-fc]?fc=o>  ^ne  matl"ices  #1  an<3  D\  can  be  computed  by  the  Levinson  algorithm,  and 
by  comparing  R\  with  (4.3)  one  finds  that 

Pjn  =  Xj>  0  <  i  <  n  and  6n  =  Xn, 

see,  e.g.,  [AG].  We  now  apply  the  recursion  formula  in  Levinson's  algorithm  backwards  in  order 
to  determine  i?i  and  D\  from  the  last  column  of  R±  and  <5ri.  Then  the  recursion  formula  is  used 
forwards  to  determine  Mq.  We  will  also  obtain  a  Hermitian  factorization  of  R  R  into  triangular 
matrices. 


Backward  Levinson  algorithm 

input:  [Pjn]?=o>^n  ;  output:  Ri,Di,  Schur  parameters  {7y}?=1  of  Mo; 
for  k  :—  n,n  —  l,n  —  2,  ...,1  do 

Ik  :=  Pofci  Pfc-i,fc-i  :=  1; 

for  j  :=  1,2,  ...,integer  part(|)  do 

solve  for  pj_iffc_i  and  /?fc-i-j,fc-i  the  linear  system  of  equations 


1       Ik 
Ik       1 


P]-l,k-l 
.Pk-l-j,k-l  . 


Pl,k 
~Pk-3,k 


&-i  :=  (&/(i  -  W))/(i  +  W); 


Levinson  recursion  for  computing  M0 

input:  Hi, JDi,{7y}y=1;  output:  {/*y}y~J; 

/i0  :=  ^o;  A*i  :=  -$o7i; 
for  A;  :=  1,2,  ...,n  -  1  do 

Mfc+i  :=  -*fc7fc+i  -  £y=i  Wy-i.fc; 


[Mj'-fcjjfcio 


Hence  Mq,R\,  and  Z?i  are  computed  in  0(n2)  arithmetic  operations  from  the  coefficients  of  x-  Let 
Rq  and  Do  denote  the  n  x  n  leading  principal  submatrices  of  R\  and  D\  respectively.  Similarly  to 
(4.4)  we  have 

RgMoRo  =  D0.  (4.5) 

1/2 

Because  Mo  is  positive  definite,  so  is  Dq.  D0      can  therefore  easily  be  computed.  We  obtain  from 
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(4.1)  -  (4.2)  and  (4.5),  with  R  :=  Dl0/2R^\ 


RHR={RT{{)T{RT1H).  (4.6) 

The  right  hand  side  of  (4.6)  is  a  Hermitian  factorization  into  triangular  matrices.  It  can  be  computed 
in  0{  n2)  arithmetic  operations  from  the  coefficients  of  X- 
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